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Abstract. The next-to-minimal supersymmetric standard model predicts the formation of 
domain walls due to the spontaneous breaking of the discrete .^-symmetry at the electroweak 
phase transition, and they collapse before the epoch of big bang nucleosynthesis if there exists 
a small bias term in the potential which explicitly breaks the discrete symmetry. Signatures 
of gravitational waves produced from these unstable domain walls are estimated and their 
parameter dependence is investigated. It is shown that the amplitude of gravitational waves 
becomes generically large in the decoupling limit, and that their frequency is low enough to 
be probed in future pulsar timing observations. 
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1 Introduction 

Supersymmetry (SUSY) is a well-motivated extension of the standard model (SM) of particle 
physics involving a solution to the hierarchy problem [see, e.g., [1, 2] for reviews]. The simplest 
supersymmetric extension of the SM is called the minimal supersymmetric standard model 
(MSSM), which consists of only SM particles along with a pair of Higgs doublets and their 
superpartners. Such a model is also motivated by the unification of the gauge forces and 
the existence of dark matter candidates. These interesting features will be investigated in 
various upcoming observations including high-energy physics, astrophysics, and cosmology 
observable. 

Aside from the rich phenomenology encountered in the MSSM, it suffers from the so- 
called /i-problem [3], which is originated from the fact that the superpotential of the MSSM 
contains a term of the form / iH u Hd (called the “^-term”) at the tree level. Here H u and H^ 
are two Higgs doublet superfields, and /r is some dimensionful parameter. Since the MSSM 
might be a low energy effective theory derived from a more fundamental theory at higher 
energies, we naively expect that the magnitude of the /i parameter is of the order of some 
cutoff scale such as the Planck scale. On the other hand, /r should be of the order of the soft 
SUSY breaking scale in order to induce the electroweak symmetry breaking appropriately. 
These considerations seem to revive the hierarchy problem in the MSSM. 

The next-to-minimal supersymmetric standard model (NMSSM) [for reviews, see [4, 
5]] is one of possible solutions to the /x-problem, which introduces an additional singlet 
superfield S and replaces the //-term with the coupling of the form A SH u Hd. Owing to this 
coupling, an appropriate magnitude of the ^ parameter is dynamically generated when the 
scalar component of the singlet superfield acquires a vacuum expectation value (VEV) at the 
electroweak phase transition. However, the introduction of the singlet superfield allows the 
existence of other dimensionful parameters at the tree level. The simplest way to prohibit 
such additional dimensionful parameters is to impose a discrete ^-symmetry, and such a 
model is called the ^-invariant NMSSM, or often abridged as the NMSSM. 

The discrete Zj, symmetry of the NMSSM is spontaneously broken when the scalar 
components of the Higgs fields as well as the singlet scalar acquire VEVs at the electroweak 
phase transition, which leads to the formation of domain walls in the universe. The existence 
of such domain walls turns out to be problematic in the standard cosmology [6], and hence it 
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was argued that the NMSSM is unfavorable from the cosmological point of view [7]. Later, 
it was pointed out that there are some possibilities to arrange the model to have a small 
explicit symmetry breaking term, which causes a late time annihilation of the walls [ 8 , 9]. 
Taking account of these possibilities, we deduce that unstable long-lived domain walls can 
exist in the early epoch of the universe. 

In this paper, we point out that a stochastic background of gravitational waves with ob¬ 
servably large amplitude can be produced from domain walls in the NMSSM. The amplitude 
of gravitational waves is determined once we specify the surface mass density of domain walls 
and their decay time [10, 11]. Here we estimate the surface mass density by solving nonlinear 
field equations for the Higgs sector in the NMSSM to find a planar domain wall solution. On 
the other hand, the peak frequency of gravitational waves solely depends on the decay time 
scale. We see that the typical frequency of the gravitational waves corresponds to the band 
which can be probed via pulsar timing observations. Another interesting finding is that the 
amplitude of gravitational waves becomes generically large in the so-called decoupling limit, 
where the couplings between the singlet field and Higgs fields become extremely small. Note 
that, due to the very weak coupling between the singlet-like states and the MSSM sector 
in the decoupling regime, the signatures such as those from the collider tend to reduce to 
those in the MSSM, and we emphasize a complementary role of the gravitational waves as a 
possible probe to distinguish the NMSSM from the MSSM. 

Signatures of gravitational waves in the NMSSM were also investigated in Ref. [12] in 
the context of the first order electroweak phase transition. Differently from that approach, in 
this paper we focus on the production of gravitational waves from domain walls, which does 
not depend on the order of the phase transition. The crucial assumption is the occurrence of 
the spontaneous breaking of a discrete symmetry followed by the late time collapse of domain 
walls, which is a definite consequence of the NMSSM with the (approximate) ^-symmetry. 

The organization of this paper is as follows. In Sec. 2, we describe the scalar potential 
in the ^-invariant NMSSM, especially focusing on its Higgs sector. In Sec. 3, we introduce 
a numerical method to analyze the structure of domain walls, and discuss their cosmological 
evolution. Then, in Sec. 4, we estimate the amplitude and frequency of gravitational waves 
produced from domain walls and discuss prospects for future observations. Finally, Sec. 5 is 
devoted to the conclusions. 

2 Higgs potential in the NMSSM 

Let us consider the supersymmetric standard model with an additional gauge singlet super¬ 
field S. As was described in Sec. 1 , we impose a discrete Z 3 symmetry, under which every 
chiral supermultiplet <I> transform as 

—>■ e 27r */ 3< F, (2.1) 

in order to forbid the presence of dimensionful parameters. The renormalizable superpotential 
of the NMSSM allowed by the Z 3 symmetry and the gauge invariance can be written as 

Hnmssm = A SH u H d + -S 3 + WVukawa: (2.2) 

where A and k are dimensionless couplings, and HYukawa corresponds to the usual Yukawa 
couplings of the quark and lepton superfields of the MSSM. In the above equations, contrac- 
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tion over gauge indices is understood. For instance, we have 

*'=(£) < 2 ' 3 > 

for SU(2)l doublet fields, which leads 

H u H d = H+H- - H° u H° d . (2.4) 


Here and hereafter we use the same symbol for the superfield and its scalar component. 
From F-term contributions of Eq. (2.2) together with D-term contributions and soft SUSY 
breaking contributions, we obtain Higgs potential at the tree level: 


U = |A (H+H- - H° u H ° d ) + kS 2 I 2 + (m 2 Hu + \X\ 2 \S\ 2 )(\H° U \ 2 + \H+\‘ 


+ K + \M \sni\H%\ 2 + \h~\ 2 ) + 


g 2 + gl 


\H° u \ 2 + \H+\ 2 -\H° d \ 2 -\H d \ 2 ) 2 


92 


+ f | H+H°* + H° u H d *\ 2 +m 2 s \S\ 2 + 


XA x (H+H d - H®H d )S + ^kA k S 3 + h.c. 


1 


(2.5) 


where g\ and g 2 are the U(l)y and SU(2)l gauge couplings, respectively, and rn 2 Hu , 

A \, and A K are dimensionful soft parameters. 

We can rotate away the component H+ by using SU (2 )l transformation, and hereafter 
we take //+ = 0. Furthermore, it is necessary to take HT = 0 in order to guarantee that the 
squared masses of charged components become positive. Therefore, in the following we only 
consider the potential for the neutral components, and simply denote them as H u = H® and 
H d = H d . The scalar potential becomes 


V=\kS 2 - XH u H d \ 2 + m 2 Hu \H u \ 2 + m\ d \H d | 2 + m 2 s \S\ 2 + |A| 2 |5| 2 (|^, a | 2 + \H d \ 2 ) 



m 2 ) 2 + 


-kA k S 3 

3 


A A x H u H d S + h.c. 


( 2 . 6 ) 


where 


2 _ 9i + 92 


2 m 2 7 

—£■- ~ 0.275 


(2.7) 


for the Z boson mass mz — 91.2GeV and the Higgs VEV v ~ 246GeV. Here we assume that 
there is no explicit CP violation, and take A, n, A x , and A K to be real. 

It is possible to take the VEV of (H u ) to be real and positive by the use of U(l)y 
transformation. In general, ( H d ) and ( S) are complex, but such an extremum becomes 
always a local maximum if there is no explicit CP violation [13]. Therefore, we can also take 
( H d ) and ( S ) to be real up to the Z 3 transformations . 1 Then, we parameterize three VEVs 
by using three real quantities: 

(H u ) = % (H d ) = % (S ) = 




yft' 


x If we parameterize general three complex VEVs as 


{Hu) = 


(Hd) = 


v/2’ 


( 5 ) = ^ 


with v u , Vd and v s being real, relevant CP-violating phases are following two combinations: 

0 = + (pd + (j>s and 5 = 3(/> s . 


( 2 . 10 ) 


(2.8) 


(2.9) 
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with v u > 0. Furthermore, we see that the potential possesses the following symmetries: A, 
k, v s —>• —A, —k, —v s and A, v d —> —A, — Vd [14]. Using these symmetries, we fix the signs of 
A and Vd to be positive. Other parameters k, v s . A x , and A K can take both signs. 

The stationary conditions of the potential dV/dH u = dV/dHd = dV/dS = 0 lead to 

Q 2 

A 2 («s + v 2 d )v u - \nv d v 2 s - Y v u( v d - v l ) + 2m % u v u - V2\A x v d v s = 0, (2.11) 

Q 2 

A 2 (u 2 + vl)v d - \kv u v 2 + y v d (v d - vl) + 2 m 2 Hd v d - V2XA x v u v s = 0 , ( 2 . 12 ) 

A 2 (vl + v 2 d )v s + 2 k 2 vI - 2 A Kv u VdV s + 2m 2 s v s - V2\A x v u v d + V2 kA k v 2 s = 0 . (2.13) 

Once we specify the values of v u , v d , and v s , the soft masses 'rri 2 Hj , m 2 Hd , and are deter¬ 
mined from Eqs. (2.ll)-(2.13) : 

^2 ^ J2 

m 2 Hu = —/r 2 - v 2 cos 2 (3 + — /i 2 cot j3 + — v 2 cos 2 (3 + fiA x cot (3, (2.14) 

u 2 A 4 

A 2 k a 2 

m 2 H = —/j 2 —— v 2 sin 2 (3 + — tan j3 - -v 2 cos 2/3 + /mA x tan/?, (2-15) 

d 2 A 4 

m 2 s = -^rv 2 - 2-^j/i 2 + ^Aku 2 sin 2/3 + A—sin2/3 — ^-/iA k , (2.16) 

Z A z 2 4 fi A 

where _ 

^V^f’ tan/3 = 5?? and V = \/ V u + v d- ( 2 - 17 ) 

Hereafter, we use the following six quantities as free parameters of the model: 

A, At, A x , A k , tan/3, and /i. (2-18) 


The definition of the /x-parameter in Eq. (2.17) implies that the value of v s becomes 
much larger than the electroweak scale for A< 1 even though /i takes a value of 0(lOO)GeV. 
Such a case is called the decoupling limit, since the coupling between particles associated 
with the singlet superfield and the MSSM sector becomes negligible for A —>• 0. In this limit, 
dominant terms in the potential (2.6) are given by 


V ~ k 2 \S\ a + m 2 s \S\ 2 + 


^ kA k S 3 + h.c. 


from which we estimate the value of v s at the minimum as 


v s 





(2.19) 


( 2 . 20 ) 


Note that the condition Af. > 9 m 2 s should be satisfied in order that the extremum given by 
Eq. (2.20) becomes deeper than the other possible minimum with v s = 0 [15]. Assuming that 

The no-go theorem of Ref. [13] ensures that these two phases take trivial values at the minimum of the 
potential. If we fix the U{1)y gauge such that v u > 0 and </>„ = 0, this fact implies that there are two more 
vacua with ((/>d,</> s ) = (2n/3,An/3) and (47r/3, 27 t/ 3) in addition to that given by Eq. (2.10). Note that it is 
also possible to choose the 17(1)v gauge such that these three degenerate vacua are represented by <j>i = 0, 
2tv/3, and 47 t/3 ( i = u,d,s), respectively. 
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the magnitudes of \A K \ and \ms\ are comparable to the soft SUSY breaking mass scale, from 
Eq. (2.20) we see that a large value of v s is compatible with k< 1. Furthermore, it can be 
shown that k and A should be comparable such that the potential does not have unrealistic 
minima [15, 16]. Therefore, the decoupling limit is given by 

A ~ k —>• 0 and v s ~ — A k /k 3> n, (2-21) 

with other dimensionful parameters being fixed. 

We emphasize that the relation [Eq. (2.21)] does not hold if v s exceeds some critical 
value, since there might exist non-renormalizable operators which dominate over k 2 |S'| 4 term 
in Eq. (2.19) when we make the value of k sufficiently small. In this regime, the value of v s is 
not given by Eq. (2.20), but determined by minimizing the scalar potential including possible 
non-renormalizable terms. Therefore, we expect that the relation [Eq. (2.21)] remains valid 
only for some intermediate range of couplings, and that the value of v s saturates to some 
finite value for a sufficiently small value of k. 


3 Domain walls in the NMSSM 


The ^-symmetry of the potential [Eq. (2.6)] is spontaneously broken when the scalar fields 
H u , Hd and S acquire VEVs at the electroweak phase transition. At that time, sheet-like 
objects called domain walls are created as a consequence of the spontaneous breaking of the 
discrete symmetry. In this section, we describe the properties of domain walls created in 
the NMSSM and some theoretical issues related to them. First, we analyze the domain wall 
solution to estimate the surface mass density in Sec. 3.1. Then, we discuss some cosmological 
aspects of NMSSM domain walls in Sec. 3.2. 


3.1 Structure of domain walls 


Since the potential is invariant under the Z% transformation [Eq. (2.1)], there exist three de¬ 
generate vacua parameterized by (( S ), ( H u ), (Hj)) = (v s /y/2, v u /y/2, Vd/V^), (v s e l27T / 3 /\/2, 
v u e l27T / 3 /V%, Vde i2n / 3 /V 2), and (v s e l4 ' K ^ 3 /y/2, v u e l4n ^ 3 /\/2, /\/2)- Domain walls are 

located around boundaries of these three vacua. 

In order to find the domain wall solution [i.e., find the spatial configuration of the scalar 
fields (S'(x), H u (x.), 77d( x )) around the core of the domain wall], we must solve nonlinear field 
equations for three complex scalar fields S, H u and Hd with the potential given by Eq. (2.6). 
It is straightforward to solve such a system numerically [7, 17]: A planar domain wall solution 
(S(z), H u (z), Hd{z)) perpendicular to the z-axis can be found from the following equations 


d 2 S _ dV d 2 H u dV d 2 H d _ dV 

~d£ ~ dS *’ dz 2 ~ OH*’ an dz 2 ~ OH*' 

with boundary conditions 


(S(z),H u (z),H d (z)) 

(S(z),H u (z),H d (z)) 


Vs Vu Vd \ 

V2’V2’V2) ’ 

•27r • 27r • 27T \ 

v s e l 3 Vue 1 3 Vd e l 3 \ 

V2 ’ V2 ’ V2 )' 


(3.1) 


(3.2) 


Here we use the globally convergent Newton method [18] to solve this boundary value prob¬ 
lem. Namely, we first specify the “initial guess” for (S(z), H u (z), H d (z)) within some finite 
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interval — L < z < L, and deform it such that Eq. (3.1) is satisfied with an appropriate 
boundary conditions specified at z = ±L. Then, we iterate this procedure until Eq. (3.1) is 
satisfied with an accuracy of 0( ICE 5 ). For an appropriate choice of the initial guess of the 
functions (S(z), H u (z), H d (z)), this algorithm converges after a few iteration steps. 

Figure 1 shows a typical configuration of the domain wall solution. From Fig. 1 (a) we 
see that the phases of scalar fields continuously change around the core of the domain walls 
centered at z = 0. We also compute the spatial distribution of the energy density for this 
solution 


Pwan(z) = 


dS_ 
dz * 


+ 


dH„ 


dz 


(*) 


+ 


dH d 


dz 


(*) 


+ V(S(z), H u (z), H d (z)) - V(v s /V2,v u /V2,v d /V2), 


(3.3) 


where we subtracted a constant term such that p wa u —>• 0 is satisfied for z —>• ±oo. As shown 
in Fig. 1 (b), the energy of the scalar fields is concentrated around the core of the domain 
wall. Finally, we estimate the surface mass density a wa n of the domain wall by integrating 
this energy density along the z-axis: 


f^waii = J dzp wa n{z). (3.4) 

For a choice of parameters used in Fig. 1, we obtain <r wa n ~ 1.9 x 10 9 GeV 3 . 




(a) (b) 

Figure 1. Spatial distribution of (a) phases and (b) energy density of scalar fields ( S(z ), H u (z), Hd{z)) 
of the planer domain wall solution for the parameters A = 0.05, n = 0.02, A\ = 150GeV, A K = 
— 150GeV, tan/? = 5 and p = 200GeV. In the left panel, <j> a , <t> u and (j>d represent the phases of S, H u 
and iJrf, respectively. Here we just show the sum of (f> u and </>d, since the relative phase between H u 
and Hd suffer from an arbitrariness due to the U(l)y gauge transformations. 

Although the precise value of cr wa u can only be calculated by using the numerical pro¬ 
cedure, we can roughly estimate the parameter dependence of <r wa n in the simple decoupling 
limit given by Eq. (2.21). Regarding the fact that the phase of the S field continuously 
changes across the core of the wall and that the potential can be approximated as Eq. (2.19) 
in the decoupling limit, we expect that the typical length scale of the spatial variation of the 
field S around the core of the wall is given by Az ~ \k,A k v s \~ 1 ^‘ 2 . Furthermore, the estima¬ 
tion shown in Eq. (2.20) implies that all three terms in the right-hand side of Eq. (2.19) are 
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comparable, and hence we roughly estimate the height of the potential around the core of 
the wall as V ~ \kA k v^\. These observations lead to the following scaling behavior of the 
surface mass density of domain walls: 

O'Wall ~ Kvl- (3.5) 

In fact, we confirmed that this behavior holds very accurately for A, k < 0(1CE 4 ) in our 
numerical results. Eq. (3.5) indicates that the surface mass density of domain walls becomes 
enormously large in the decoupling limit. We will discuss the observational consequence of 
this property in Sec. 4. 


3.2 Domain wall problem and its solution 


Let us turn our attention to the cosmological evolution of domain walls. After the formation 
of domain walls, two kinds of forces act on them. One is the tension force, which smooths 
out small scale structures to straighten the wall, and the other is the friction force caused 
by the interactions between particles in the thermal plasma and Higgs fields constituting 
the wall. The friction force is efficiently induced by particles whose masses are comparable 
to the electroweak scale, since their interactions with the Higgs fields are relatively strong. 
However, the number density of these particles drops exponentially as the temperature of the 
plasma T decreases, and eventually the friction force becomes negligible. The friction effect 
for the NMSSM domain walls was quantitatively estimated in Ref. [7], and it was shown that 
this effect becomes irrelevant for T < 0(0.1 l)GeV. 

Once the friction force becomes negligible, the dynamics of domain walls is mostly 
determined by the tension force, which acts like px ~ oVaii/i^waii for a patch of the walls 
whose typical curvature radius is given by R W aii- In this case, it is known that the evolution 
of domain walls is described by the so-called scaling solution, in which the typical length 
scales of the system are given by the Hubble radius. In particular, the energy density of 
domain walls in the scaling regime is given by [ 10 ] 


Pwall — A 


Cwall 

t ’ 


(3.6) 


where A is a dimensionless coefficient called the area parameter. By comparing a naive 
estimation p wa n ~ cr W aii/^waii with Eq. (3.6), we estimate the typical curvature radius of 
the walls in the scaling regime as R wa ii ~ t/A, and hence the tension force is given by 
PT ~ Acr wa \\/t. 

In the literature, the scaling behavior of domain walls [Eq. (3.6)] is confirmed both 
numerically and analytically [10, 19]. The result of the numerical simulation for the model of 
a real scalar field with Z 2 symmetry indicates A — 0.8 ± 0.1 [ 10 ], while other simulations for 
the model with N degenerate vacua show that the value of A increases proportionally with 
N [20]. Since we are interested in the model with three degenerate vacua (i.e. Z 3 invariant 
NMSSM), we simply estimate the value of A as A ~ 0.8 x (3/2) = 1.2. 

If the domain walls evolve according to the scaling solution, the decrement of the energy 
density /5 wa u oc 1/t is slower than that of dusts oc a(t )~ 3 and radiations oc a(t)~ 4 , where a(t) 
is the scale factor of the universe. This fact implies that domain walls are likely to dominate 
the energy density of the universe at later times. Such an existence of domain walls leads to 
many problems in the standard cosmological scenario [ 6 ]. 

A simple solution to the domain wall problem is to introduce an additional term AV 
which explicitly breaks the discrete symmetry in the potential [ 6 , 21]. Let us call this addi¬ 
tional term a bias. We note that the magnitude of the bias term should be much less than 
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that of the potential around the core of domain walls (AV <C V) such that the discrete 
./^-symmetry holds approximately. If such a bias term exists in the potential, it breaks the 
degeneracy between different vacua, and the energy difference between neighboring domains 
causes the volume pressure py ~ AV acting on the walls. Domain walls start to collapse 
when this pressure py exceeds the tension force px ~ Aa %va _\\/t : and we can estimate their 
decay time td ec from the condition py ~ px- 


^dec — 


“^^"wall 

AV ' 


(3.7) 


If the decay of domain walls occurs before they overdose the universe, we have p wa ii{tdec) < 
Pent (idee) = 3H 2 (fdec)/87rG, where /9crit(£dec) is the critical density of the universe at the time 
tdeo H(tdec) is the Hubble parameter at that time, and G is the Newton’s constant. From 
this requirement, we obtain 


cT wa ii < 2.93 x 10 4 TeV 3 A -1 


/ O.lsec 
\ ^dec 


(3.8) 


Even if domain walls decay before they overdose the universe, the consideration of big 
bang nucleosynthesis (BBN) leads to another constraint if they are sufficiently long-lived. 
For domain walls existing around the epoch of BBN, the ratio between their energy density 
and the entropy density s is given by 


~2x 10 -7 GeV A 
s 


OVall \ 
ITeV 3 / 



(3.9) 


It is expected that particles radiated from the walls are likely to destroy light elements created 
during the epoch of BBN. For the abundance of radiated particles, which is comparable 
with Eq. (3.9), we require fd ec 1$ O.Olsec in order to avoid the constraints on hadronic 
energy injection in the standard BBN scenario [22], This condition gives a lower limit on the 
magnitude of the bias term: 


AV > 6.6 x 10" 2 MeV 4 A wal1 , . (3.10) 

Several remarks on the origin of the bias term [Eq. (3.10)] were made in the literature. 
A naive expectation is that the discrete symmetry is slightly broken due to gravity, and 
that a Planck-suppressed higher dimensional operator is responsible for the bias term [23]. 
However, it turned out that this possibility is unfavorable in the NMSSM [7], since such a 
Planck-suppressed interaction radiatively induces a tadpole operator [24], which revives the 
hierarchy problem because of the large VEV of the singlet field. A solution to this problem 
was proposed in Ref. [8]: It was pointed out that a small tadpole term of the form 

AP - £m 3 oft S + h.c. (3.11) 


can be generated by imposing various additional symmetries to constrain the form of non- 
renormalizable interactions. Here, m so f t is the soft SUSY breaking mass, and the dimension¬ 
less coefficient £ corresponds to the loop suppression factor. The term given by Eq. (3.11) is 
small enough not to destabilize the singlet field, and it acts as a bias term to remove domain 
walls. 2 

2 We note that there are other various attempts to avoid tadpole and domain wall problem [25-27]. These 

models also introduce additional symmetries to forbid a large tadpole term. Since the form of the Lagrangian 

densities of these models becomes different from that of the Z 3 -invariant NMSSM [Eq. (2.2)], which does not 
lead to the formation of domain walls, we do not consider these possibilities in this paper. 
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In addition to the solution described above, there is another possibility, where a coupling 
between the singlet field and additional vector-like matters which are charged under the SM 
SU(3)c gauge group or a hidden QCD gauge group is introduced [9]. In the presence of 
such a coupling, the ^-symmetry becomes anomalous for the QCD, and the bias term which 
destabilizes the domain walls is generated due to the instanton effect [28]. In this case, the 
magnitude of the bias term is estimated as 

AV ~ A 4 , (3.12) 

where A corresponds to the QCD scale Aqcd ~ 0(lOO)MeV if vector-like exotics are charged 
under SU(3)c, or a scale at which the hidden gauge force becomes strong if they are charged 
under the hidden gauge group. 

Once the magnitude of the bias term is specified as Eq. (3.11) or Eq. (3.12), we can 
estimate the decay time of domain walls according to Eq. (3.7). Since the origin of the bias 
term is highly model dependent, we take Q ec as a f ree parameter in the rest of this paper. 


4 Gravitational waves from domain walls 


From the discussion in the previous section, we see that domain walls exist after the elec- 
troweak phase transition, and that they collapse before the epoch of BBN f^ec O.Olsec 
due to the effect of the bias term. In this section, we estimate gravitational wave signatures 
produced from these unstable domain walls [10, 11, 29]. 

The cosmological evolution of domain wall networks is accompanied by various violent 
processes of collisions and annihilations of them. During these processes, a part of the energy 
stored in them is released as gravitational waves, which leads to a stochastic gravitational 
wave background observed today. The amplitude of gravitational waves at the cosmic time 
t for a certain frequency / can be characterized by the following quantity 


O m — 1 dp gw (t) 

gwl j Pcrit(i) dlnf 


(4.1) 


where p c rit(Q and p gw (f) are the critical density of the universe at the time t and the energy 
density of gravitational waves at that time, respectively. 

The spectrum of gravitational waves from domain walls is estimated in Refs. [10, 11] 
by the use of 3-dimensional lattice simulation of the scalar field. According to the recent 
detailed study [10], the spectrum of gravitational waves from domain walls has a peak at the 
frequency corresponding to the Hubble radius at the time idee- The peak amplitude of the 
gravitational waves at that time is given by [10] 3 


(Qec)peak 


STregwG 2 ^ 2 ^ 

3ET 2 (tdec) ’ 


(4.2) 


where e gw is an efficiency parameter of the gravitational radiation. The result of the numerical 
simulation gives e gw ~ 0.7T0.4 [10]. Assuming that the gravitational radiation is terminated 

3 In addition to the contribution from the dynamics of domain walls, there exists a secondary background 
due to the effort of the scalar fields to maintain the scaling defect [30]. The spectrum of such a secondary 
background turns out to be scale invariant, extending over a wider frequency range. The amplitude of the 
scale-invariant background figw can be estimated from the VEV of the scalar held, and in the NMSSM case we 
obtain Bgw(tdec) ~ (u s /Afpi) 4 , where Mp ~ 1.22x 10 19 GeV is the Planck mass. Since this amplitude is smaller 
than Eq. (4.2) by a factor of ~ (H(td ec ) / A K ) 2 , hereafter we do not consider this secondary contribution. 
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in the radiation dominated era, we estimate the amplitude of the gravitational waves at the 
present time to as 

^ gw h 2 (t 0 ) ^n R h 2 (j-j Ogw(tdec) 

/in 75 \ 1/73 

~ 3.45 x lO" 5 x l Qgw(t dec ), (4.3) 

where ^IrH 2 = 4.15 x 10~ 5 is the density parameter of radiations at the present time, h is the 
present value of the Hubble parameter in the unit of 100km-sec _1 Mpc _1 , gso = 3.91 is the 
effective degrees of freedom for the entropy density at the present time, and go = 3.36 and 
g * are the number of radiation degrees of freedom at the present time and t^ ec , respectively. 
Substituting Eqs. (3.7) and (4.2) into Eq. (4.3), we obtain 

ttgw^oWk - 5 - 20 X 10 " 20 X e ~gw-4 4 (lTeV 3 ) { ^AV ) ' ^ 


The peak frequency is given by the Hubble parameter at the decay time: 


f(t oW = 

a(t 0 ) 

~ 3.99 x 10 _9 Hz A~ 1/2 


ITeV 3 \ 1/2 
°wall J 


/ AH \ 1/2 

V IMeV 4 / 


(4.5) 


We see that 12 gw /i 2 (fo) pe ak and /(to)peak are determined by two parameters: <7 wa u and 
AV [or those appearing in the explicit from of the bias term, i.e., Eq. (3.11) or Eq. (3.12)]. 
Furthermore, AV is related to the decay time of domain walls via Eq. (3.7), and hence 
Eqs. (4.4) and (4.5) are rephrased as 


^gw/i 2 (to)peak^l.20x 10 17 xe gw .4 4 ( ^ ) (lTeV 3 ) (o.otsec) 
/(io)peak - 1-02 X 10 _9 Hz 


'0.01sec\ x / 2 
^dec / 


(4.6) 

(4.7) 


As was noted in Sec. 3.2, here we use t^ ec as a free parameter instead of AH. 

In order to estimate the peak amplitude for a given set of the NMSSM parameters 
[Eq. (2.18)], we must evaluate the value of the surface mass density of domain walls <j wa n 
appearing in Eq. (4.6). Here, we use the numerical method introduced in Sec. 3.1 to estimate 
<7 W aii• Figure 2 shows the estimated peak amplitude in two-dimensional parameter space of 
(A, k). In this figure, we also plotted some constraints on the tree-level Higgs potential: One is 
the condition that the mass matrix of seven Higgs scalars (three CP-even, two CP-odd, and a 
pair of charged states) do not have any negative eigenvalues, and the other is the requirement 
that the Higgs potential does not lead to unrealistic minima. We used the analytic expression 
for unrealistic minima obtained in Ref. [16], and checked whether the depth of the desired 
minimum with ((£), (H u ), (Hd)) = (v s /y/2, v u /y/2, Vd/V%) is deeper than that of possible 
unrealistic minima for given values of the parameters (A, k). 4 

4 An additional constraint might be imposed in order to avoid charged and colored vacua [15]. We do not 

put this constraint explicitly in Fig. 2, since it also depends on other MSSM parameters such as soft masses 

for sleptons and squarks. 
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Figure 2. The peak amplitude of gravitational waves estimated from Eq. (4.6) with the assumption 
of tdec = O.Olsec in the (A, k) plane. Blue solid lines correspond to contours for R gw ^ 2 (to)peak = 10 -9 , 
10 -12 , and 10 -15 . Here we use <?* = 10.75, A = 1.2 and e gw = 0.7 for numerical coefficients appearing 
in Eq. (4.6). Other theoretical parameters are fixed as A\ = 150GeV, A K = —150GeV, tan/3 = 5 
and n = 200GeV. The dark green region corresponds to the parameters at which the mass matrix of 
Higgs scalars has some negative eigenvalues, and the light green region corresponds to the parameters 
at which some unrealistic minimum becomes deeper than the desired minimum. The pink region 
corresponds to the parameter space where the amplitude of gravitational waves exceeds the present 
limit fig w h 2 < 10 -8 from PTA observations. The gray region corresponds to the parameter space 
where domain walls dominate the energy density of the universe before tdec [he. Eq. (3.8)]. 


The increase in the peak amplitude of gravitational waves for A <C 1 is straightforwardly 
understood from the relation in the decoupling limit [Eq. (3.5)]. This relation implies that 
O gw ^ 2 (io)peak oc K 2 v & s t 2 Aec oc K 2 A~ 6 ^ 6 ^dec i n the limit A —» 0.° We see that the decoupling limit 
is relevant to the observations, since it generically predicts a large amplitude of gravitational 
waves. Regarding this fact, in Table 1 we choose some bench mark points in the parameter 
space, whose predictions can be tested in future experimental research. 

There are various ongoing and planned gravitational wave experiments, which aim to de¬ 
tect stochastic gravitational wave backgrounds as well as astrophysical sources. Among them, 
ground-based interferometers such as (Advanced) LIGO [32], VIRGO [33], KAGRA [34], and 
ET [35] will probe frequency range around / ~ 100Hz. Space-borne interferometers such as 
eLISA [36] and DECIGO [37] are also planned to probe lower frequency ranges around O.lrnHz 
to 10Hz. In addition to these experiments, gravitational wave backgrounds can be probed by 
using the Pulsar Timing Array (PTA) [38], which is sensitive to the frequency range around 
1CT 9 -10 -8 Hz. Recently, three projects including PPTA, EPTA, and NANOGrav reported 
upper limit on the stochastic gravitational wave background of fl gw /i 2 < 0(1) x ICE 8 [39]. 
As shown in Fig. 2, this present upper limit already excludes some parameter region for 

s Note that the coupling A cannot be arbitrary small with keeping the relation in the decoupling limit 
[Eq. (2.21)], as was discussed at the end of Sec. 2. This fact implies that we cannot extract arbitrarily large 
amounts of the energy of gravitational waves by simply keeping reducing the coupling A. 
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Table 1 . Three bench mark points and estimated surface mass density of domain walls, peak am¬ 
plitude of gravitational waves, and its frequency. Here we used g* = 10.75 for = 10~ 2 sec and 
g * = 68.8 for idee = 10 _6 sec [31] to estimate n gw /i 2 (io) P eak- 



I 

II 

III 

A 

5 x nr 4 

5 x nr 3 

5 x 10 " 9 

K 

2 x 1 (T 4 

2 x 10~ 3 

2 x 10 -6 

A\ 

150GeV 

150GeV 

150GeV 

a k 

— 150GeV 

— 150GeV 

— 150GeV 

tan (3 

5 

5 

5 


200GeV 

200GeV 

200GeV 

Mec 

10 ” 2 sec 

10 2 sec 

10 6 sec 

<2"wall 

1.96 x 10 4 TeV 3 

1.96 x 10 2 TeV 3 

1.96 x 10 s TeV 3 

Qgwh (to)peak 

4.66 x 10 -9 

4.66 x 10 “ 13 

2.51 x 10 -9 

/ (i'O)peak 

1.02 x 10 9 Hz 

1.02 x 10~ 9 Hz 

1.02 x 10“ 7 Hz 


tdec = O.Olsec. The future PTA observations such as FAST [40] and SKA [41] will further 
improve the sensitivities. 

Figure 3 shows the schematics of the spectrum of gravitational waves predicted from 
three points enumerated in Table 1 and sensitivities of planned experiments. The recent 
result of numerical simulations [ 10 ] indicates that the spectrum of gravitational waves from 
domain walls scales like S4 gw oc / 3 for / < / pea k and fi gw oc / _1 for / > f pea k, and we 
also plot this expectation in the figure. We see that the point I leads to the peak at the 
frequency / ~ 10~ 9 Hz, which can be probed in PTA observations such as SKA. On the other 
hand, the prediction of the point II indicates that the amplitude gets smaller as the values 
of A and k increase. Interestingly, the spectrum of gravitational waves predicted from the 
point III might be probed by DECIGO as well as SKA. From these illustrations, we expect 
that the future PTA observations and space-borne interferometers can probe the signature 
of gravitational waves from NMSSM domain walls. 

5 Conclusions 

In this paper, we estimated the gravitational wave signatures from domain walls created at the 
electroweak phase transition in the NMSSM. The surface mass density of domain walls < 7 wa u is 
estimated by solving nonlinear equations for Higgs fields numerically. By using the estimated 
value of cT wail, the amplitude of relic gravitational waves is calculated via Eq. (4.6) and its 
parameter dependence is explored. We find that a large amplitude of gravitational waves, 
which scales as Al g w^p Ca k oc is predicted in the decoupling limit. Furthermore, the 

peak frequency is determined by the decay time of domain walls [Eq. (4.7)], which is related 
to the magnitude of the bias term AV via Eq. (3.7). Such gravitational wave signatures can 
be probed in the future experiments with improved sensitivities. 

We note that the decay time of domain walls f<jec is restricted to a certain interval, 
from the epoch of the electroweak phase transition to that of BBN. The typical frequency of 
gravitational waves corresponding to such time scales is very low, which is suitable for PTA 
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ET 

eLISA 

DECIGO. 

SKA . 

point I - 

point II -- 

point III - 


Figure 3. Sensitivity curves of planned detectors (dotted lines) and the spectrum of the gravitational 
waves produced from domain walls in the NMSSM (solid lines) for various choices of parameters (points 
I, II, and III) given in Tabe 1. The sensitivity curve for ET is based on [42], and those for eLISA and 
(Ultimate) DECIGO are based on [43, 44] [see also Ref. [10]]. The sensitivity of SKA is taken from 
Ref. [45], 


observations. Since the amplitude of gravitational waves becomes large in the decoupling 
limit, we can use the results of PTA observations as a lower limit on the couplings (A and 
k) appearing in the Higgs sector of the NMSSM. For instance, the present PTA bound 
Og W h 2 < 0(1O- 8 ) already puts a lower limit A > 10 -4 -10 -3 for t <± ec = O.Olsec (see Fig. 2), 
while this constraint becomes irrelevant if the decay time of domain walls t^ec is earlier than 
0(O.Ol)sec. Future observations will probe a wider range of the parameters in the weak 
coupling regime, which can be used to complement the results of other experimental studies 
of the NMSSM. 
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